######### Figures #####

rm(list=ls())
dev.off()

#Libraries
library(data.table)
library(lfe)
library(multiwayvcov)
library(lmtest)

library(readxl)
library(foreign)

library(dplyr)
library(raster)
library(stringr)
library(abjutils)

library(stargazer)

source("paysdetat_functions.R")

##### DATA ########

load("livingstandards_data.RData")
load("rebellion_data.RData")
load("grievance_data.RData")

## Remove pays d'imposition and ambiguous province Venaissin ##
dat <- reb[generalite!="Venaissin",]
reb <- reb[imposition==0,]

te <- te[imposition==0,]
n <- n[imposition==0,]
p <- p[imposition==0,]

controls <- c("wheat","rugged","distparis","distriver","distsea","distforeignborder","cities_10_50km_1400","latitude","longitude","latitude*longitude")

###### Figure 4: Balance checks ###########

pdf("balancechecks.pdf",height=9,width=6)

vars <- c("wheat","rugged","distparis","distriver","distsea","distborder","distcity10k_1400","cities_10_50km_1400","distintendance")
labvars <- c("wheat","rugged","dist. Paris","dist. river","dist. sea","dist. border","dist city","N cities","dist. intendance")

par(mfrow=c(5,1))
par(mar=c(2,12,1,3))

xs <- seq(10,50,1)

for(i in 1:length(vars)){
  coefs <- lapply(xs*1000, dist="distestate",FUN=rdmodel, datause=population[imposition==0,], loc="nearET + x + y + x*y", cont="x",depvar=vars[i])
  coefs <- do.call(rbind,coefs)
  lim <- 1.5*max(coefs[,1]+1.96*coefs[,2],abs(coefs[,1]-1.96*coefs[,2]))
  plot(xs,coefs[,1], ylim=c(-lim,lim), ylab="",yaxt='n')
  axis(4)
  axis(2,at=0,label=labvars[i],las=1,tick=F,cex.axis=1.5)
  segments(xs,coefs[,1]+1.96*coefs[,2],xs,coefs[,1]-1.96*coefs[,2])
  abline(h=0,lty=2)
}

dev.off()


###### Figure 5: RD: Local rule and dissatisfaction with royal taxation ######

pdf("rd_tax.pdf", width=9, height=5)

par(mfrow=c(2,1))
par(mar=c(2,9,3,3))

vars <- c("taxall_s")
labvars <- c("% grievances\nabout taxation\n(Third Estate)","% grievances\nabout taxation\n(Nobility)")

xs <- 20:50
rdplot(dat=te,var="taxall_s",labvar="% grievances\nabout taxation\n(Third Estate)",cont=controls, ylim=c(-2,2))
rdplot(dat=n,var="taxall_s",labvar="% grievances\nabout taxation\n(Nobility)",cont=controls, ylim=c(-2,2))

vars <- c("sharefiscal","sharefiscalnosalt")
labvars <- c("% fiscal\nrebellions","% fiscal\nrebellions\nexcluding\nsalt tax")

xs <- 10:50
rdplot(dat=reb,var="sharefiscal",labvar="% fiscal\nrebellions",cont=controls)
rdplot(dat=reb,var="sharefiscalnosalt",labvar="% fiscal\nrebellions\nexcluding\nsalt tax",cont=controls)

dev.off()


###### Figure 6: RD: grievances about the seigneurial regime #########

pdf("rd_elites.pdf",width=7,height=3)

par(mfrow=c(2,1))
par(mar=c(2,10,3,3))
xs <- seq(20,50,1)
rdplot(dat=te,var="seign_s",labvar="% grievances\nabout seign. regime",cont=controls, ylim=c(-2,2))
rdplot(dat=te[nearET!="S",],var="seign_s",labvar="% grievances\nabout seign. regime\n(weak Third Estate)",cont=controls, ylim=c(-2,2))

xs <- seq(10,50,1)
rdplot(reb,"shareseign",labvar="% rebellions\nagainst seign system")
rdplot(reb,"sharelocal",labvar="% rebellions\nagainst local elites")

rdplot(reb[nearET!="S",],"shareseign",labvar="% rebellions\nagainst seign system\n(weak Third Estate)")
rdplot(reb[nearET!="S",],"sharelocal",labvar="% rebellions\nagainst local elites\n(weak Third Estate)")

dev.off()



###### Figure B2: Characteristics of general lists corresponding to the village sample #####
te_rebin20km <- te_rebin20km[imposition==0,]
te_rebin20km[,treat:=ifelse(b_id %in% p[imposition==0,b_id],1,0)]

vars <- c("wheat","rugged","cities_10_50km_1400","distriver","distsea","distforeignborder","distintendance","maxpop1793","n_griev","rebelsum")
labs <- c("wheat","ruggedness","urban","dist. river","dist. sea","dist. land border","dist. intendance","max pop 1793","N Third Estate grievances","N rebellions")
cbind(vars,labs)

tab <- lapply(vars, FUN=function(x){felm(as.formula(paste("scale(",x,") ~ treat | 0 | 0 | generalite",sep="")), data=te_rebin20km) %>% summary %>% .$coefficients %>% .["treat",]}) %>% 
  do.call(what=rbind,.) %>% as.data.table

coefs <- tab$Estimate
ses <- tab$`Cluster s.e.`

pdf("balance_villages.pdf")

par(mfrow=c(1,1))
par(mar=c(5,12,3,3))

xs <- 1:length(coefs)
plot(coefs,xs, xlim=c(-1.5,1.5), yaxt='n', ylab="")
segments(coefs-1.96*ses,xs,coefs+1.96*ses,xs)
abline(v=0, lty=2)
axis(2,at=xs,labels=labs,las=1)

dev.off()




###### Figure C.1: Relationship between 1789 Third Estate grievances and 1700-1788 popular rebellions #########

f1 <- felm(taxall_s ~ sharefiscal+ncoms+scale(log(maxpop1793))+log(1+pop1793) + wheat + distparis + cities_10_50km_1400 | 0 | 0 | generalite, data=te)
f2 <- felm(taxall_s ~ sharefiscal+ncoms +scale(log(maxpop1793))+log(1+pop1793)+wheat+distparis + cities_10_50km_1400 | 0 | 0 | generalite, data=d2[maxpop1793<7871,])
f3 <- felm(taxall_s ~ sharefiscal+ncoms +scale(log(maxpop1793))+log(1+pop1793)+wheat+distparis + cities_10_50km_1400 | 0 | 0 | generalite, data=d2[maxpop1793>=7871,])

f4 <- felm(seignsse_s ~ sharelocal+ncoms +scale(log(maxpop1793))+log(1+pop1793)+wheat+distparis + cities_10_50km_1400 | 0 | 0 | generalite, data=d2)
f5 <- felm(seignsse_s ~ sharelocal+ncoms +scale(log(maxpop1793))+log(1+pop1793)+wheat+distparis + cities_10_50km_1400 | 0 | 0 | generalite, data=d2[maxpop1793<7871,]) 
f6 <- felm(seignsse_s ~ sharelocal+ncoms +scale(log(maxpop1793))+log(1+pop1793)+wheat+distparis + cities_10_50km_1400 | 0 | 0 | generalite, data=d2[maxpop1793>=7871,]) 

pdf("griev_rebel_cor.pdf",width=9, height=6)

models <- list(f1,f2,f3,f4,f5,f6)

coefs <- lapply(models,function(x){x$coefficients[2]}) %>% unlist
ses <- lapply(models,function(x){x$cse[2]}) %>% unlist
xs <- c(2,3,4,6,7,8)

plot(xs, coefs, ylim=c(-.5,.7),xlim=c(1,9),xaxt='n',xlab="",pch=21,bg=rep(c("gray","white","black"),2))
segments(xs,coefs-1.96*ses,xs,coefs+1.96*ses)
abline(h=0, lty=2)
legend('topright',title="1789 districts:",c("all","rural","urban"),pch=21, pt.bg=c("gray","white","black"))
text(1,-.2,labels=c("Dep. var: % grievances about taxation"),adj=0)
text(1,-.3,labels=c("Indep. var: % fiscal rebellions"),adj=0)

text(5.5,-.2,labels=c("Dep. var: % grievances about seign. regime"),adj=0)
text(5.5,-.3,labels=c("Indep. var: % rebellions against local elites"),adj=0)

abline(v=5)

dev.off()

###### Figure D.1. Proportion of peasant grievances in urban Third Estate and Nobility lists ######
load("grievlevel_data.RData")

#How similar are village lists to general TE lists and Nobility lists ?
similarity <- function(d1,d2,pat){
  
  b_id <- unique(intersect(d1$doc_id_b_no,d2$doc_id_b_no))
  
  out <- rep(NA,length(b_id))
  
  for(i in 1:length(b_id)){
    #Select grievances for each district
    x1 <- d1$griev_c_f4[d1$doc_id_b_no==b_id[i]] #parish
    x2 <- d2$griev_c_f4[d2$doc_id_b_no==b_id[i]] #bailliage
    x1 <- grep(pat,x1,value=T) 
    x2 <- grep(pat,x2,value=T)
    
    inter <- length(unique(intersect(x1,x2))) #how many common
    #union <- length(unique(union(x1,x2)))  #how many
    village <- length(unique(x1))  #how many
    #Denominator: N unique grievances by both in baillage
    #But why is this "unique"? should not it be all grievances?
    out[i] <- inter/village
  }
  return(out)
}

d2 <- t_bygriev[t_bygriev$doc_id_b_no %in% te$doc_id_b_no,]
vt <- similarity(p_bygriev[imposition==0,],d2 ,"")
vt_tax <- similarity(p_bygriev[imposition==0,],d2,"^GTA")
vt_seign <- similarity(p_bygriev[imposition==0,],d2,"^SSE")

vn <- similarity(p_bygriev[imposition==0,],n_bygriev,"")
vn_tax <- similarity(p_bygriev[imposition==0,],n_bygriev,"^GTA")
vn_seign <- similarity(p_bygriev[imposition==0,],n_bygriev,"^SSE")

pdf("boxplot_villages_otherlists.pdf")

boxplot(cbind(vt,vt_tax,vt_seign,vn,vn_tax,vn_seign), names=rep(c("all","tax","seign."),2), col=c(0,0,0,'lightblue','lightblue','lightblue'))
legend('topright',c("general Third Estate","Nobility"), pt.bg=c(0,'lightblue'),pch=22)

dev.off()


###### Figure E3: RD: impact of local rule on per capita land tax in 1802 ######

pdf("rd_landtaxpc1802.pdf", height=4, width=9)

controls <- c("wheat","rugged","distparis","distriver","distsea","distforeignborder","distintendance","cities_10_50km_1400")

vars <- c("contfonc1802pc")
labvars <- c("1802 arrondissement\nper capita land tax")

par(mfrow=c(1,1))
par(mar=c(2,14,1,3))

xs <- 20:50
coefs <- lapply(xs*1000, dist="distestate",FUN=rdmodel, datause=cont[imposition==0,], loc="nearET + x + y + x*y", cont=controls,depvar=vars)
coefs <- do.call(rbind,coefs)
lim <- 1.5*max(coefs[,1]+1.96*coefs[,2],abs(coefs[,1]-1.96*coefs[,2]))
plot(xs,coefs[,1], ylim=c(-lim,lim), ylab="",yaxt='n')
axis(4)
axis(2,at=0,label=labvars,las=1,tick=F,cex.axis=1.5)
segments(xs,coefs[,1]+1.96*coefs[,2],xs,coefs[,1]-1.96*coefs[,2])
abline(h=0,lty=2)

dev.off()

###### Figure E4: OLS and RD: grievances about various aspects of taxation ########

pdf("rd_taxtype.pdf", height=9)

vars <- c("taxall_s","indtax2_s","dirtax2_s","taille_s","univtax_s","taxadm_s","gtadire_s","taxadvpriv_s")
labvars <- c("taxation","indirect","direct","taille","universal","tax collection","tax distribution","tax privilege")

par(mfrow=c(length(vars)/2,1))
par(mar=c(2,10,3,3))

xs <- c(20:50)
for(i in 1:length(vars)){
  rdplot(te,var=vars[i],labvar=labvars[i], ylim=c(-2,2))
}
for(i in 1:length(vars)){
  rdplot(n,var=vars[i],labvar=labvars[i], ylim=c(-2,2))
}

dev.off()



###### Figure F1: balance checks (OLS) #######
#Geographic / pop variables
vars <- c("wheat","rugged","distparis","distriver","distsea","distborder","distcity10k_1400","cities_10_50km_1400","distintendance")
coefs <- list()
for(i in 1:length(vars)){
  coefs[[i]] <- lapply(10^10, dist="distestate",FUN=rdmodel, datause=population[imposition==0,], loc="nearET + coordx + coordy + coordx*coordy", cont="coordx",depvar=vars[i]) %>% unlist
}
coefs <- do.call(rbind, coefs)

pdf("balanceols.pdf")

labvars <- c("wheat","rugged","dist. Paris","dist. river","dist. sea","dist. border","dist city","N cities","dist. intendance")
cbind(vars,labvars)

par(mfrow=c(1,1))
par(mar=c(3,10,2,2))
xs <- rev(1:nrow(coefs))
plot(coefs[,1],xs, xlim=c(-.75,.75), ylim=c(0,nrow(coefs)+1), yaxt='n', ylab="")
abline(v=0, lty=2)
segments(coefs[,1]+1.96*coefs[,2],xs,coefs[,1]-1.96*coefs[,2],xs)
axis(2,at=xs,labels=labvars, las=2)

dev.off()


###### Figure F2: Robustness checks: drop border sections (rebellion data) #####

pdf("rd_conflict_dropborders.pdf",height=9,width=7)

par(mfrow=c(5,1))
par(mar=c(2,10,3,2))

xs <- 10:50
rdplot(reb[nearET!="BRETAGNE",], "sharefiscal","drop Brittany",main="Fiscal rebellions")
rdplot(reb[nearET!="N",], "sharefiscal","drop Northern",main="")
rdplot(reb[nearET!="BOURGOGNE",], "sharefiscal","drop Burgundy",main="")
rdplot(reb[nearET!="S",], "sharefiscal","drop Southern",main="")
rdplot(reb[nearET!="SW",], "sharefiscal","drop Pyrenean",main="")

rdplot(reb[nearET!="BRETAGNE",], "sharelocal","drop Brittany",main="Local elites")
rdplot(reb[nearET!="N",], "sharelocal","drop Northern",main="")
rdplot(reb[nearET!="BOURGOGNE",], "sharelocal","drop Burgundy",main="")
rdplot(reb[nearET!="S",], "sharelocal","drop Southern",main="")
rdplot(reb[nearET!="SW",], "sharelocal","drop Pyrenean",main="")

dev.off()


###### Figure F3: Robustness checks: drop border sections (grievance data) #####

pdf("rd_griev_dropborders.pdf",height=9,width=7)

controls <- c("wheat","rugged","distparis","distriver","distsea","distforeignborder","distintendance","cities_10_50km_1400")

vars <- c("taxall_s","seign_s")
labvars <- vars
nearETs <- c("BRETAGNE","NORD","BOURGOGNE","S","SW")
labs <- c("drop Brittany","drop Northern","drop Burgundy","drop Southern","drop Pyrenean")

par(mfrow=c(length(nearETs),1))
par(mar=c(2,9,3,3))

xs <- c(20:50)
for(i in 1:length(nearETs)){
  rdplot(te[nearET!=nearETs[i],], var="taxall_s", labvar=labs[i], cont=controls, ylim=c(-2.5,2.5))
}

xs <- c(20:50)
for(i in 1:length(nearETs)){
  rdplot(te[nearET!=nearETs[i],], var="seign_s", labvar=labs[i], cont=controls, ylim=c(-2.5,2.5))
}

dev.off()


###### Figure F4: Robustness checks: control for type of primary source rebellion data #########

clust <- "generalite"
sourcevars <- c("archnat","archdep","municip","bib")

coefs <- list()

depvars1 <- c("fiscalrebel","fiscalnosaltrebel","localrebel","seignrebel")
coefs[[1]] <- lapply(depvars1,FUN=rdmodel,
                     datause=reb,treat="until1789",clust="generalite",cont=c("pop1793",sourcevars,controls), summary=F) %>%
  lapply(., FUN=function(x){x["until1789",]}) %>% do.call(what="rbind") %>% as.data.frame

depvars2 <- c("fiscalpop","fiscalnosaltpop","localpop","seignpop")
coefs[[2]] <- lapply(depvars2,FUN=rdmodel,
                     datause=reb,treat="until1789",clust="generalite",cont=c(sourcevars,controls), summary=F) %>% 
  lapply(., FUN=function(x){x["until1789",]}) %>% do.call(what="rbind") %>% as.data.frame
coefs[[3]] <- lapply(depvars1, FUN=rdmodel, datause=reb, cont=c("pop1793",controls,sourcevars), loc="latitude + longitude + latitude*longitude + nearET", bw=30*1000) %>% do.call(rbind,.) %>% as.data.frame
coefs[[4]] <- lapply(depvars2, FUN=rdmodel, datause=reb, cont=c(controls,sourcevars), loc="latitude + longitude + latitude*longitude + nearET", bw=30*1000) %>% do.call(rbind,.) %>% as.data.frame

names(coefs[[3]]) <- names(coefs[[1]])
names(coefs[[4]]) <- names(coefs[[1]])

coefs <- do.call(rbind,coefs) %>% as.data.table

coefs$depvar <- rep(c(depvars1,depvars2),2)
coefs$ols <- c(rep(1,8),rep(0,8))
coefs$typevar <- gsub("rebel|pop","",coefs$depvar)

coefs <- coefs[order(typevar),]

pdf("coefplot_rebellion_sources.pdf", height=5)

par(mar=c(5,3,3,2))

xs <- 1:nrow(coefs)

plot(xs,coefs$Estimate, ylim=c(-1,1), pch=c(1,2,3,4),  xaxt='n',xlab="",ylab="Impact of local rule")
segments(xs,coefs$Estimate+1.96*coefs$`Std. Error`,xs,coefs$Estimate-1.96*coefs$`Std. Error`)

abline(h=seq(-1,1,.2), col='gray', lty=2)
points(xs,coefs$Estimate, ylim=c(-1,1), pch=c(1,2,3,4))

abline(h=0, lty=2)
abline(v=c(4.5,8.5,12.5))

legend("topright",c("OLS: % rebellion (1)","OLS: rebellion p.c. (2)","RD 30 km: % rebellion (3)","RD 30 km: rebellion p.c. (4)"),pch=1:4,bg="white")

mtext("Fiscal",side=1,at=2.5)
mtext("Fiscal\nw/o smuggling",side=1,at=6.5, line=1)
mtext("Local\nelites",side=1,at=10.5, line=1)
mtext("Seign.\nelites",side=1,at=14.5,line=1)

dev.off()


###### Figure F5: RD: local rule and grievances (count outcome) ######

pdf("rd_griev_count_dum.pdf",height=7,width=7)

par(mar=c(2,12,3,3))
par(mfrow=c(2,1))
xs <- seq(20,50,1)
rdplot(te,"taxall","N grievances\nabout taxation\n(Third Estate)",main="", cont=c(controls,"n_griev"), ylim=c(-2.5,2.5))
rdplot(n,"taxall","N grievances\nabout taxation\n(Nobility)",main="", cont=c(controls,"n_griev"),  ylim=c(-2.5,2.5))
rdplot(te,"seignsse","N grievances\nabout seign. regime",main="", cont=c(controls,"n_griev"),  ylim=c(-2.5,2.5))
rdplot(te[nearET!="S",],"seignsse","N grievances\nabout seign. regime\n(weak Third Estate)",main="", cont=c(controls,"n_griev"),  ylim=c(-2.5,2.5))

dev.off()

###### Figure F6: RD: local rule and fiscal rebellions (count and binary outcome)########

pdf("rd_conflict_tax_count_dum.pdf",height=9,width=7)
xs <- 10:50
par(mfrow=c(4,1))
par(mar=c(2,12,3,2))
rdplot(reb, "fiscal","Count fiscal rebellion",main="")
rdplot(reb[city_code!=44109,], "fiscal","Count fiscal rebellion\nExclude Nantes",main="")
rdplot(reb, "fiscalnosalt","Count fiscal rebellion\nNo smuggling",main="")
rdplot(reb[city_code!=44109,], "fiscalnosalt","Count fiscal rebellion\nNo smuggling\nExclude Nantes",main="")
rdplot(reb, "fiscal1","Dummy fiscal rebellion",main="")
rdplot(reb, "fiscalnosalt1","Dummy fiscal rebellion\nNo smuggling",main="")
rdplot(reb, "fiscalpop","Fiscal rebellion\nper capita",main="")
rdplot(reb, "fiscalnosaltpop","Fiscal rebellion\nper capita\nNo smuggling",main="")
dev.off()



###### Figure F7: RD: local rule and rebellions against local elites (count and binary outcome) ########

pdf("rd_conflict_elites_count_dum.pdf",height=9,width=7)

xs <- 10:50
par(mfrow=c(3,1))
par(mar=c(2,12,3,2))
rdplot(reb, "seign","Count rebel\nagainst seign. elites",main="")
rdplot(reb, "seign1","Dummy rebel\nagainst seign. elites",main="")
rdplot(reb, "seignpop","Rebel against\nseign. elites\nper capita",main="")
rdplot(reb, "local","Count rebel against\nlocal elites count",main="")
rdplot(reb, "local1","Dummy rebel against\nlocal elites",main="")
rdplot(reb, "localpop","Rebel against local elites\n per capita",main="")

dev.off()

